Torsten Hothorn, Frank Bretz and Peter Westfall (2008), Simultaneous Inference in General Parametric Models. Biometrical Journal, 50(3), 346-363; See vignette("generalsiminf", package = "multcomp").
人工データ
set.seed(300)
count <- c(rpois(30, 10) + floor(rnorm(30, 0, 3)), rpois(30, 10) +
floor(rnorm(30, 0, 3)), rpois(30, 14) + floor(rnorm(30, 0, 3)))
treat <- factor(c(rep("A", 30), rep("B", 30), rep("C", 30)))
d <- data.frame(count, treat)
head(d)
## count treat
## 1 14 A
## 2 14 A
## 3 14 A
## 4 7 A
## 5 7 A
## 6 18 A
ポアソン回帰
res <- glm(count ~ treat, d, family = poisson)
summary(res)
##
## Call:
## glm(formula = count ~ treat, family = poisson, data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.090 -0.923 -0.032 0.940 2.991
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2756 0.0585 38.88 <2e-16 ***
## treatB 0.0370 0.0820 0.45 0.65
## treatC 0.3847 0.0759 5.07 4e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 207.47 on 89 degrees of freedom
## Residual deviance: 174.69 on 87 degrees of freedom
## AIC: 555.9
##
## Number of Fisher Scoring iterations: 4
##
Tukeyの方法で多重比較
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: splines
res.t <- glht(res, linfct = mcp(treat = "Tukey"))
summary(res.t)
##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: glm(formula = count ~ treat, family = poisson, data = d)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## B - A == 0 0.0370 0.0820 0.45 0.89
## C - A == 0 0.3847 0.0759 5.07 <1e-05 ***
## C - B == 0 0.3477 0.0750 4.63 1e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
同時信頼区間による多重比較
res.ci <- confint(res.t)
res.ci$confint <- exp(res.ci$confint)
plot(res.ci, xlab = "relative risks")
abline(v = 1, lty = "dashed")